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ABSTRACT 

The solidification of a dilute alloy (bismuth-tin) under Bridgman crystal growth conditions is 
investigated. Computations are performed in two dimensions with a uniform grid. The simulation 
includes the species concentration, temperature and flow fields, as well as conduction in the 
ampoule. Fully transient simulations have been performed, with no simplifying steady state 
approximations. Results are obtained under microgravity conditions for pure bismuth, and for Bi-0. 1 
at.%Sn and Bi-1.0 at.%Sn alloys, and compared with experimental results obtained from crystals 
grown in the microgravity environment of space. For the Bi-1 .0at.%Sn case the results indicate that a 
secondary convective cell, driven by solutal gradients, forms near the interface. The magnitude of the 
velocities in this cell increases with time, causing increasing solute segregation at the solid/liquid 
interface. The concentration-dependence of the melting temperature is incorporated in the model for 
the Bi-1.0 at.%Sn alloy. Satisfactory correspondence is obtained between the predicted and 
experimental results in terms of solute concentrations in the solidified crystal. 
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NOMENCLATURE 

c p specific heat at constant pressure 

C species concentration 

D species diffusion coefficient 

f volume fraction 

g mass fraction 

g acceleration due to gravity 

Gr Grashof number, gJ3-r(TH-Tc)H 3 /v 2 

Gr s solutal Grashof number, gB c CoH 3 /v 2 

h ampoule thickness (outside radius - inside radius) 

ho reference enthalpy = ? H/p 

H ampoule diameter; reference length 

k thermal conductivity 

k p segregation coefficient 

L length of simulation domain 

La translating zone length 

Le Lewis number, a/D 

Pr Prandtl number, v/a 

t time 

T temperature 

u, v velocities in x and y directions 

Greek Symbols 

a thermal diffusivity 

B expansion coefficient 

? T temperature difference, Th - Tc 

? H enthalpy of freezing 

v kinematic viscosity 

? radial segregation, (C m ax ■ C m j n )/Caverage 

? density 

Subscripts 

0 initial condition 
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C cold furnace temperature condition 
H hot furnace temperature condition 
L liquid 

m at solidification front 
S solid 
w ampoule wall 
Superscripts 
~ vector 

* reference quantity 

INTRODUCTION 

The synthesis of advanced materials, especially for electronics and biomedical applications, 
demands high-quality crystals. The compositional uniformity (and hence the quality) of such crystals 
can be profoundly influenced by the transport phenomena which occur in the melt during 
solidification. The primary transport mechanism causing these deleterious effects is natural 
convection. The low-gravity environment of space offers an opportunity to suppress the strength of 
this natural convection. Hence there is a great deal of interest in the study of directional 
solidification of crystals in space. 

The MEPHISTO project (Abbaschian et al., 1992) is a collaborative program of space 
experiments aimed at understanding the fundamental processes involved in crystal growth. The 
space-borne experimental apparatus is a Bridgman-type furnace with an isothermal hot zone, an 
isothermal chill zone, and an insulated gradient zone. The furnace contains three ingots of Bi-Sn or 
Sn-Bi binary alloy inside fused silica ampoules with a maximum 6 mm inner diameter and a 10 mm 
outer diameter. All three samples are solidified simultaneously under identical thermal conditions. 
After flight, the samples are extracted and analyzed. Four MEPHISTO space experiments have 
taken place; the most recent, MEPHISTO-4, flew in November 1997. The MEPHISTO-2 and -4 
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experiments examined the faceted solidification of bismuth doped with tin (Bi-Sn). 

The experimental data from the MEPHISTO-2 mission experiments (with a Bi-0. 1 at.% Sn alloy) 
resulted in a greater understanding of the dominant role of interface kinetics on morphological 
stability . Stability phenomena were observed that had not been previously predicted by theory or 
measured in terrestrial experiments (Abbaschian, 1 996). The MEPHISTO-4 experiments built on to 
the previous results, and used a Bi-1 .0 at.% Sn alloy. 

The MEPHISTO project includes a program of computational modeling of the crystal growth 
process. In particular, the role of convection, which is crucial to a complete understanding of the 
process, is to be investigated. Since accurate experimental determination of convection in metallic 
melts is very difficult to achieve, due to the opacity and chemical reactivity of the melts, convective 
levels are determined numerically. Furthermore, the computational models themselves are to be 
improved by a process involving prediction of, and comparison with, the experimental results. The 
aim of this procedure is to develop effective fully transient computer simulations of fluid flow 
related effects. Previous computations of Bridgman growth had been limited to steady-state growth 
in succinonitrile, a widely used transparent phase change material with properties analogous to 
metallic materials (Yao and de Groh, 1993; de Groh and Yao, 1994). 

Convection effects at microgravity levels were modeled by means of a transient, 2D FIDAP 
finite-element model by Yao et al. (1997). A fixed-grid approach was adopted, with the enthalpy 
method being employed to model the phase change. Temporal averaging was used for the apparent 
heat capacity in the discretized equations. Due to computational difficulties introduced by the small 
partition coefficient for Bi-Sn, the presence of solute was ignored in Yao et al. (1997). Preliminary 
scaling arguments by de Groh and Nelson (1994) implied that solutal convection effects on solute 
segregation may be significant. However, only recently has it been possible to include solutal 
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convection into numerical simulations involving phase change for MEPHISTO-4. This was due to 
the difficulties concerning convergence with front-tracking methods as well as those imposed by the 
low partition coefficient for Bi-Sn (Yao et al., 1995; Yao et al., 1997). 

Many simulations of Bridgman crystal growth processes, both under terrestrial and microgravity 
conditions, are available in the literature. The majority of these simulations can be classified as 
pseudo steady state models. The key assumption in such models is that a “steady-state” mode of 
alloy solidification exists, i.e. the concentration of the dopant in the solid which forms at the 
interface is equal to the initial dopant concentration in the liquid (Kurz and Fisher, 1989). Such 
models vary in complexity from simple 2D analyses that consider the interface to be flat (Alexander 
et al., 1989; Yao et al., 1995; Simpson et al., 1998a) to much more complex formulations that are 
able to handle interface curvature and wall conduction (Adornato and Brown, 1987) and fully 3D 
simulations (Liang and Lan, 1996). However, because of the low partition coefficient for Bi-Sn 
alloys, a steady-state mode of solidification is never achieved during the MEPHISTO experiments. 
Thus pseudo steady state models are not appropriate; recently developed fully transient simulations, 
such as those of Simpson and Garimella (2000) and Timchenko et al. (2000) need to be employed 
and improved in order to faithfully model this process. 

The computational modeling presented in this paper is intended to examine the effects of 
thermosolutal natural convection on the MEPHISTO-4 solidification experiments. This will be 
achieved using a fully transient 2D model, which includes most of the effects of binary alloy 
solidification, convection driven by both thermal and solutal gradients, distinct thermal properties in 
the solid and liquid phases and the effects of interface curvature. Compared to Simpson et al. 
(1998b), the solution scheme used in this work is improved in terms of the following aspects. 
Conduction through the ampoule wall is considered. The dependence of melting temperature on 
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concentration is included for the richer alloy. It must be stressed that for the growth rates and 
applied thermal gradients encountered during the highly controlled directional solidification 
experiments modeled in this paper, the liquid/solid interface was experimentally found to be stable 
and non-dendritic like that for a pure metal. Consequently, the solid/liquid interface modeled in our 
numerical analysis is also stable and non-dendritic. 

MATHEMATICAL FORMULATION 

The problem under consideration is the directional solidification of a binary alloy by the 
Bridgman process, as shown schematically in Fig. 1. The gravity vector is perpendicular to the 
furnace axis (horizontal Bridgman growth configuration). The melt region is considered to be a 
viscous Newtonian fluid subject to thermosolutal convection. Thermophysical properties are 
considered as constant but distinct for the solid and liquid phases. Density variations are considered 
to be subject to the Boussinesq approximation. The equation for conservation of momentum and 
mass are 

^ + (yxu)xu=-p T {T-T c )g-p c (C-C 0 )g-^ + xK 2 u ( 1 ) 

Po 

V-M=0 (2) 

For calculating nondimensional constants (such as thermal Rayleigh number, Ra) the ampoule 
inside diameter H is selected to be the reference length. The characteristic time and velocity become 
t (= H 2 /a) and v* (= H/t = a/H). 

The governing equation for the conservation of energy is 

a{t ' Cl) aT + ^«^(“ r) =v-(tw)+5(r,c t ) (3) 

It will be seen below that effective heat capacity A and source term B may be specified in this way 
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for calculating phase-change. An initial temperature equal to Th is applied throughout the flow field. 
The temperatures at the x = 0 and x = L walls are set to be T H and Tq, respectively. 

The thermal boundary conditions along y = • (H/2 + h) are a function of time and are shown 
schematically in Fig. 1 . There is a translating zone (considered anAadiabatic® zone if the temperature 
profile is unknown) between the hot and cold regions of the furnace, in which the temperature 
linearly increases from the cold furnace temperature to the hot furnace temperature. The melting 
temperature of the material occurs somewhere within this zone, which translates with time at a 
constant x-velocity, known as the translation velocity, u,. This is what facilitates the directional 
growth of the crystal. Defining the x location where the translating zone meets the cold furnace 
temperature zone to be at xa(0, the boundary condition for temperature may be expressed as 

T c , x<x A (t ) 

y = t = \t c + at^^-, * x (/)<;x<(x,(r )+l a ) (4) 

1 J L A 

T h , ( x A (t) + L A )<x 

In principle, the solution of the energy equation (3) coupled with the solution of the momentum and 
mass conservation equations (1 and 2) would yield the temperature and velocity distribution 
throughout the simulation domain. However, the problem of modeling the physics of the 
propagation of the solidification front and determining its location remains to be addressed. We 
choose to do this by employing the phase-transformation model of Zeng and Faghri (1994). In this 
model, apparent heat capacity A and source term B are given as 
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A(T,C L )= pc p + a{T,C L )^r 

o T 

B(T,C l ) = -V ■ ( ph 0 u)+ V • [pf s (h 0 + (c pL 


-^)T)u s ]-a(T,C L ) ^ 

dC L dt 


a{T,C L )= p 


c ps )T "*■ 


( Pi ~ Ps)^jT_ + PlK 
Pl ~ f L (P L ~ Ps') 


c p = c p lSl + Cpsgs 
P ~ Psfs + P L ft 


(5a) 


For the present study , the density of each phase is assumed to be equal, and there is no dispersion of 
solid phase moving in the liquid. For the results shown where the melting temperature is considered 
dependent on concentration, the full expressions for A and B (Eq. 5a) simplify to become: 

A(T,C,) = pc„+a(T)^ 

Ol 

£(r ’ c ' )= - a(r >l^ (5b) 

a{T) = p[(c pL - CpS )T + h Q ] 

For the pure bismuth and dilute alloy results, concentration-dependence on temperature is neglected, 
and the expressions for A and B simplify further to 


A(T)^pc r+ a(X)^ 

oT 

B(T ) = 0 

a(7-) = p[(< v -c^)r+* 0 | 

The equation for conservation of solute throughout the computational domain is 
^+V-(«Q) =DV 2 C l +S(T,C l ) 


(5c) 


( 6 ) 


This equation is analogous to the energy equation. We impose an initial solute concentration 
throughout the solution domain. At the boundaries no solute may exit the solution domain. Thus, 
t-0 Cl = Co 

x = 0, L • C[/»x = 0 / 7 n 
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y = 0,H 


*Ci/»y = 0 


Again, in principle, solution of equation (6) along with energy equation (4) and the fluid 
velocities all subject to the relevant boundary and initial conditions is enough to determine the 
solute, temperature and velocity values throughout the solution domain. However, the more general 
problem involving phase change demands that solute redistribution at the advancing solid/liquid 
interface be addressed. At the interface, new solid material forms at a concentration lower than at 
the adjacent liquid, in accordance with the partition coefficient (Kurz and Fisher, 1989): 

= ^p^L ( 8 ) 

For the Bi-Sn system under consideration in this work, the partition coefficient k p has a value of 
0.029. This small value indicates that the solid forms at a concentration very poor in solute, and thus 
there is rapid and significant build-up of solute in the liquid region ahead of the interface. Source 
term S accounts for solute rejection at the interface into the bulk fluid resulting from the effect of the 
solid material forming at the lower concentration. Following the work of Swaminathan and Voller 
(1997) and Voller et al. (1989) source term S may be written as the differential 

S(T,C L ) = '* f £ L) +k p C L %S- (9) 

NUMERICAL ANALYSIS 

The numerical scheme is explained in detail in Simpson and Garimella (1 998, 2000); only salient 
details are provided here. The computational domain is primarily discretized using regularly spaced 
finite difference mesh points. Superimposed on this grid are finite volumes which are used for the 
solution of the energy and species concentration equations. The finite volume centers are staggered 
with respect to the finite difference mesh point locations. Vorticity , velocity and vector potential are 
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calculated at the finite difference mesh points. Temperature and solute concentration are evaluated 
using the finite volumes. 

For solving velocities, the vorticity-streamfiinction representation of the Navier-Stokes and 
continuity equations (1 and 2) is used. The numerical approach applied here is modified from a 
program (Timchenko et ah, 1997) written for the solution of natural convection in a rectangular 
cavity. The essential details are that the discretized equation for the vorticity transport equation is 
solved using an Alternating Direction Implicit scheme (Samarskii and Andreyev, 1963). The 
discretized equation for streamfunction is solved using the conjugate gradient method. Once the 
values of streamfunction are known, the nodal velocities can be determined. Boundary conditions 
are explicitly applied on the advancing solid/liquid interface which is arbitrarily oriented and so the 
boundary conditions require special handling. The location and slope of the interface is calculated 
using a Hirt and Nichols (1981) type front reconstruction; the boundary conditions may then be 
applied once the slope and position of the front are known. 

The energy (4) and solute conservation (6) equations are discretized using the finite volume 
mesh, with an upwind scheme incorporated for the treatment of convective heat fluxes. The 
discretized equations are solved using iterative solvers - either Gauss-Seidel iteration or Stone’s 
strongly implicit procedure (Ferziger and Peric, 1996). For the general case where melting 
temperature is dependent on interface concentration, the concentration and energy equations are 
highly coupled and are solved simultaneously. 

RESULTS AND DISCUSSION 

Validation 

A grid-independence study was undertaken to determine the appropriate spatial and temporal 


10 


discretization scheme to be used. This study involved examining the solution fields from a test 
matrix of simulations performed using four different discretization schemes. Four mesh sizes of 
increasing spatial and temporal refinement were examined. The mesh selected on the basis of this 
grid independence study was a regularly spaced 300 • 35 mesh, with 25 of the 35 y-direction cells 
being in the sample and the remaining 1 0 in the ampoule walls. The time step size was ? t = 0. 1 338 
s. Detailed results from the grid-independence study, including plots of key field variables, are 
available in Simpson and Garimella (2000). The key points are: 

• The thermal field results were almost identical for all four meshes. 

• The difference between the velocities predicted by the mesh adopted in this study and those 
from a finer mesh was less than 3%. 

• Particular attention was given to the concentration values in the solidified material. The 
concentrations in the solid are completely dependent on the solutal, thermal and flow fields 
in the melt, and so, are very sensitive to any changes in these fields. Concentration traces in 
the solid for the mesh employed in this study had a maximum discrepancy of only 3.2% (with 
the majority of values being less than 1%) compared to results at the finer mesh. 

These results make it clear that (a) all the field variables converge to a single result as the mesh 
becomes increasingly fine, and (b) the mesh employed in this work is fine enough to obtain 
meaningful results while keeping computational costs reasonable. 

Pure Bismuth 

Simulations for the Bridgman crystal growth of pure bismuth were performed. The cold and hot 
furnace temperatures were Tc = 50* C and Th = 700* C respectively. Thermophysical properties 
from Yao et al. (1995) were used, and are shown in Table 1. Key nondimensional parameters 
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describing thermal transport are Gr = 1 13.6 and Pr = 0.01 144. For this pure Bismuth case, there is 
no solute, so the dilute model for phase chance (Eq. 5c) was used. 

In order to start each of the crystal growth simulations, the following procedure was carried out. 
The initial position of the translating zone was flush with the x = 0 wall. This zone is immobilized 
for the first 3,000 time steps. During this time, the velocity and concentration field solution schemes 
are switched off while solid rapidly chills in the portion of the translating zone which is lower than 
the melting temperature. This new state is then taken to be at time t = 0. After this, simulations 
proceed with the entire solution scheme enabled and the insulated zone moving at the translation 
velocity u t . 

Figure 2(a) shows a plot of velocity vectors and isotherms after 3,000 s have elapsed. Isotherms 
are shown in the ampoule wall as well as in the liquid and solid bismuth. The thick line at x • 16 
mm represents the solid/liquid interface. The dominant feature of this plot is the counter-clockwise 
convective cell in the translating zone (1 6 < x < 40 mm). For clarity, velocity vectors are plotted on 
every third mesh point in the x-direction in this figure. The maximum velocity in this convective 
cell is 3.47 pm/s at (24.75, • 1 .8 mm). The ratio of this velocity to the translation velocity, u t , is 1 .03, 
which compares to the value of 1 . 1 2 found from the finite-element, variable-property simulation of 
the same process by Yao et al. (1997). Velocities in the negative y-direction at the interface are 
constrained and concentrated by the presence of the solid whereas velocities in the positive y- 
direction at the opposite (hot) end of the translating zone are much more diffuse. 

The other main feature of this plot is the isotherms throughout the solution domain. On the 
outside edge of the ampoule, the (imposed) linear temperature profile within the translating zone is 
witnessed by the regularly spaced isotherms on the outer edge. In the low-conductivity ampoule the 
isotherms are dramatically distorted, and the thermal field on the inside of the ampoule wall differs 
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from that imposed on the outside. For x • 24 mm, the temperature in the ampoule is greater than the 
applied temperature (the isotherms distort sharply to the left); for x • 24 mm, the converse is true. In 
the bismuth region the isotherms exhibit the same trend, and have a gentle, crescent-shaped 
curvature. The isotherms and the interface appear to be symmetric about the centerline (y = 0) and 
thus have not been influenced by convective transport in the melt. The solid-liquid interface, which 
is an isotherm, is discemibly curved such that the solid is concave. The total deflection of the 
interface is about 0.5 mm which compares to a value of -0.46 mm found in Yao et al. (1997). This 
curvature is primarily a result of the difference in thermophysical properties for the solid and liquid 
phases. Translation of the thermal boundary condition also contributes slightly to this effect. 

Figures 2(b) and 2(c) are plots of velocity vectors and isotherms at 6,000 and 9,000 s, 
respectively. The progression of the translating zone (and thus of the solidification front) is obvious 
from these figures. After 6,000 s, the front is at x • 26 mm, and the nature of the convective motions 
is unchanged. The magnitude of the maximum velocities at these later times is identical to that at the 
previous time (Fig. 2a). This indicates that end effects are not influencing the results. The curvature 
of the interface also remains the same. 

Bi» 0.1 aU% Sn alloy (MEPHISTO-2) 

The growth of a Bi # 0.1 at.% Sn (0.08185 vol.% Sn) alloy was considered next. The velocity, 
thermal and solute fields were solved subject to the dilute alloy model (Eq. 5c), and the 
dimensionless parameters for solute transport are Gr s = 0.349 and Le - 3074. Simulations were 
performed using the same spatial discretization and time step as for pure bismuth. 

The nature of the convective flow for this alloy (results not shown) is similar to that for pure 
bismuth. A single counter-clockwise thermally driven convective cell dominates the domain. The 
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maximum u- velocity is slightly less than for pure bismuth, and is at the same location. However, the 
magnitude of the maximum v-velocity is • 1.77 pm/s at (17.75, 0 mm) acting in the negative y 
direction, which is lower in magnitude than for pure bismuth (• 1.87 pm/s). Solute rejected at the 
interface acts to oppose the thermally driven convective motion, resulting in retarded velocities near 
the interface. The thermal field is not distorted by the action of convection. As time progresses to 
6,000 and 9,000 s, the maximum v-velocity decreases to • 1 .72 pm/s and • 1 .65 pm/s respectively, as 
a result of continuing solute rejection at the interface. 

Traces of solute concentration across the height of the solidified alloy in the domain are shown in 
Fig. 3. Three different x-locations (16.18, 26.30, 36.41 mm) are considered, corresponding to the 
location of the interface at y = 0 at times of 3,000, 6,000 and 9,000 s, respectively. To aid in 
visualizing the influence of convection, results for a simulation involving diffusion only (no 
convection) are superimposed on this plot. Note that these pure diffusion results indicate that 
significant radial segregation arises due to interface curvature (see Adomato and Brown, 1 987; Liang 
and Lan, 1 996) such that the solute concentrations at the centerline are larger than those at the edges. 
The level of curvature-induced radial segregation remains almost constant, with values of ? = 29. 1 % 
(based on the mean concentration) at x = 16.18 mm and 29.4% at x = 26.30 and 36.41 mm. The 
inclusion of convection into the calculations causes the following effects. Solute levels are increased 
for low values of y and decreased for higher values of y. The y-location for the maximum solute 
concentration is also shifted from the centerline to y = • 1.23 mm for all the traces. This is a 
consequence of the convective flow sweeping solute Adown« the interface in the direction of 
decreasing y and away from the interface in the direction of increasing x. The maximum values for 
solute concentration are C max /C 0 = 0.30, 0.50, 0.65 for the traces at x = 16.18, 26.30 and 36.41 mm, 
respectively. The radial segregation values increase slightly from ? = 50.1% at x = 16.18 mm to 
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56.3% at x = 36.41 mm. 


Traces of solute concentrations in the solid and melt at the different times are shown in Fig. 4. 
These traces run longitudinally along the ampoule for three different heights of y= • 2.88, 0 and 2.88 
mm. Note that the interface location on the centerline (y = 0) trace lags slightly behind the value for 
the other traces since the interface is curved. Solute build-up near the interface, and the exponential 
profile characteristic of binary-alloy solute rejection (Smith et al., 1955) are clearly evident. The 
concentration at the centerline increases from C/Co = 9.62 at 3,000 s to 20.79 at 9,000 s. The impact 
of convection can also be seen in this plot. Warm, solute-poor fluid impinges on the top section of 
the interface. This fluid, cooled by the interface, falls toward the bottom wall and then returns to the 
bulk. This convective transport causes a thinner solute boundary layer and lower interface 
concentration value for y = 2.88 mm and, conversely, a thicker solute boundary layer and higher 
interface concentration value for y = • 2.88 mm. The interface concentrations along the centerline 
remain larger than those found at the periphery for all the times shown. This effect is due to 
interface curvature (see Fig. 3). The axial segregation in the solidified crystal is seen to increase with 
time; the radial segregation evident in Fig. 3 is also seen in Fig. 4. 

For the low-concentration alloy, the solute concentration profiles are governed by solute rejection 
at the (curved) interface and by thermal convection only. The action of solutal convection is limited 
at these low concentrations. The velocity field indicates that thermal convection decreases slightly 
with time. 

Bi • 1.0 at.% Sn alloy (MEPHISTO-4) 

A richer alloy was simulated next. The initial liquid composition was taken to be at a uniform 
1.0 at.% Sn (0.8185 vol.% Sn) increasing the value of Gr s to 3.49. The dilute model (Eq. 5c) was 
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retained for this initial investigation; the validity of this is commented upon later. 

Velocity vectors and isotherms at three different times are shown in Fig. 5. At 3,000 s (Fig. 5a) 
the velocities indicate a primary counter-clockwise convective cell in the translating zone, along with 
a weak secondary, clockwise cell driven by solute gradients adjacent to the interface. The two-cell 
convective motion is in contrast to that observed for pure bismuth and for the more dilute alloy. The 
maximum v-velocity in the secondary convective cell near the interface at (16.5, 0 mm) is 0.57 pm/s. 
The nature of the maximum velocities for this case is also different due to the higher levels of solute 
rejection and the presence of the secondary cell. The maximum u- velocity for this case is slightly 
lower than that for pure bismuth. The maximum v-velocity is identical to the value for pure bismuth 
but is located further from the interface (x= 19.25 versus x= 17.75 mm for the pure and Bi-0.1 at.% 
Sn cases). 

At a later time (Fig. 5b), the secondary convective cell has increased in size and strength as solute 
- the driving force for this cell - continues to build up at the interface. The maximum v-velocity in 
the secondary convective cell near the interface at (27.00, 0.48 mm) is 1.75 pm/s. Later in the 
growth process (Fig. 5c), the front has advanced to x • 36 mm, and the secondary convective cell has 

become quite strong. The maximum v-velocity in the secondary cell (and in the domain) is 2.88 
pm/s at (37.00, 0.36 mm). 

Traces of solute concentration across the height of the solidified alloy, along with pure diffusion 
results, are shown in Fig. 6. This plot provides a stark contrast with the results shown for the more 
dilute alloy in Fig. 3; the formation and growth of the secondary cell causes segregation to occur in 
the opposite sense (such that concentration values near the top of the domain are higher than at the 
bottom). At x = 1 6. 1 8 mm, the concentrations are close to those for pure diffusion; the maximum 
value is C/C 0 = 0.28 at y = 0.059 mm with a segregation value of ? = 28.9%. This serendipitous 
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result is a consequence of segregation arising from the influence of the secondary cell balancing the 
segregation caused by the thermally driven cell when solidification has proceeded to this location. 
At lower values of x (t < 3,000 s), segregation is such that the maximum value occurs at y < 0 
(similar to the Bi-0.1 at.% Sn alloy result, see Fig. 3) while at higher values of x (t > 3,000 s) this 
maximum is shifted to a more positive y value. The traces at x = 26.30 and 36.41 mm show that 
increasing solutal convection leads to increasing radial segregation. The maximum values are C/Co 
= 0.51 and 0.75 with radial segregation levels of ? = 54.7 and 82.2%, respectively. This behavior is 
different from the case shown in Fig. 3, where the convective field remained almost steady with 
respect to the interface resulting in only a slight increase in segregation. 

Figure 7 is a plot of longitudinal solute concentrations at domain heights of y =• 2.88, 0 and 2.88 
mm for the Bi-1 .0 at.% Sn case at the three different times. As for the case shown in Fig. 4 there is 
significant solute buildup, with liquid solute concentrations C/Co at the interface along the centerline 
being 9.70, 15.59 and 19.20, respectively. Note, however, that this corresponds to much higher 
concentrations than for the case shown in Fig. 4 due to the higher initial concentration Co- Near the 
interface, the solute concentrations in the liquid are higher at y = 2.88 mm than at y = • 2.88 mm. 
Further from the interface, but still within the solute boundary layer, concentration increases with 
decreasing y. This is due to the interaction of the two convective cells. Close to the interface, the 
secondary cell circulates solute-rich material within the boundary layer, acting to accumulate solute- 
rich material in this region. Further from the interface (see Fig. 5), the primary convective cell 
influences the solute distribution, sweeping solute from the bottom of the domain into the bulk and 
thickening the solute boundary layer in this region. In general, the solute boundary layers are more 
compact than for the dilute alloy case at each corresponding time. Note that concentration values of 
this magnitude would be sufficient to lower the melting temperature of the alloy by a significant 
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amount (Simpson et al., 1998a). This implies that the constant melting temperature assumption 
made in these simulations is not valid. This effect will be addressed in the next subsection. 

For the higher alloy concentrations, solutal convection plays a much larger role, as can be seen in 
the concentration profiles discussed above. As growth proceeds, the level of solutal convection 
grows, with a corresponding increase in radial segregation. The maximum concentration increases 
much more rapidly than the minimum value. 

Figure 8 is a plot of average solute concentrations in the solid at 3,000, 6,000 and 9,000 s for 
both the dilute and rich alloy. Averaging is performed across the height of the solidified material (y 
= • 3 to 3 mm). For comparison purposes, average concentrations for the pure diffusion case as well 
as those calculated from an analytical expression for solid solute concentrations under ID, pure 
diffusion growth conditions (Smith et al., 1955) are also shown. The pure diffusion results exhibit 
good agreement with the 1 D analytical result; the largest difference between these results is 1 .2% at 

3.000 s. The lowest average solid concentrations are found for the dilute alloy (Bi-0.1 at.%Sn). 
After 3,000 s the dilute alloy concentration is 1.3% less than for pure diffusion. With time, this 
deficit increases to 2.2% after 6,000 s and 2. 1 % after 9,000 s. The average solute concentrations for 
the richer alloy (Bi-1.0 at.% Sn) lie between those found for the pure diffusion and Bi-0.1 at.% Sn 
results. These concentrations are 0.76, 1 .2 and 1 .7% lower than the pure diffusion values after 3,000, 

6.000 and 9,000 s, respectively. For the dilute alloy, the single convective cell acts, on average, to 
increase the size of the solute boundary layer (Fig. 5) by transporting solute from the interface out 
into the melt (Fig. 3). This results in a corresponding decrease in the average solid concentration 
levels when compared to the pure diffusion case, for a given y-location. For the richer alloy, the 
solute boundary layers are shaped by a more complicated convection pattern (Fig. 6). In general, the 
boundary layers are more compact than for the dilute alloy case, resulting in slightly higher average 
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concentrations. The difference between the average concentration found for the pure diffusion and 
Bi-1.0 at.%Sn cases increases with increasing solutal convection levels. 

Comparison with Experiment 

An actual solidification process carried out during the MEPHISTO-4 space experiment was 
modeled using the rich-alloy model (Eq. 5b) which accounts for the effect of concentration- 
dependent melting temperature. The main difference in this simulation compared to those above is 
that the solidification cycle selected was conducted at a furnace translation rate of 1 1 .997 (im/s for a 
duration of 334 seconds, resulting in a total furnace translation of 4 mm. The gravity level for this 
case was reduced to 1 |ig, which is more typical of residual gravity levels during the experiment (de 
Groh and Nelson, 1 994). Under these conditions, the process is diffusion-dominated although there 
is some discernible convection-induced segregation (Simpson and Garimella, 2000). 

Figure 9 is a plot of concentration values along the centerline of the sample. The square symbols 
are experimental data points obtained via a post-flight microprobe analysis of the solidified sample 
(REZA/HENRY: PLEASE PROVIDE A REFERENCE FOR THE POST-FLIGHT DATA 
REPORT). Although there is scatter in the experimental data, a clear trend of increasing 
concentration values can be discerned. The solidification process depicted starts at a location 143.5 
mm along the ampoule; on Fig. 9, this location is taken as the origin x = 0 at t = 0. Then, in order to 
grow crystal from the melt, the furnaces are moved at a constant velocity. As the furnaces begin to 
move and new crystal is formed, solute begins to pile up at the interface. As the concentration at the 
interface increases, the local melting temperature decreases. Hence the front slowly “melts back” 
with respect to the furnaces and thus moves at a speed that is slower than the furnace translation 
velocity. After 334 seconds have elapsed, the furnaces have traversed 4 mm, however, the 
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sohd/hquid interface has moved only 3.1 mm, and so only 3.1 mm of ciystal is grown. The results 
from the numerical simulation, using the rich alloy model (Eq. 5b) that includes the effect of 
concentration-dependent melting temperature, are superimposed on this plot as the thick solid line. 
The experimental and numerical results are shown for a transverse location of y = 0. These results 
agree reasonably well with the experimental data both in terms of the predicted values of 
concentration and in terms of predicting the total amount of material grown. Results of a simulation 
using the dilute model (Eq 5c, no concentration-dependence of melting temperature) are also shown 
as a dashed line on the figure. This simulation clearly cannot handle the “melt back” of the interface 

and is therefore not able to predict the correct final location of the solid-liquid interface or accurate 
values of concentration. 

The effect of the inclusion of concentration-dependent melting temperature in the calculations is 
to predict less net solute rejection into the liquid at a given time, since the front moves at an overall 
slower rate. If this concentration-dependence was incorporated into the results considered in Figs. 5, 

6 and 7 for Bi-1.0 at.% Sn above, the primaiy differences would be twofold. First, less net solute 
would be rejected, so that the gradients in the liquid region near the interface would be less steep, 
resulting in lower levels of convection and hence lower levels of segregation at a given time. The 
other mam effect would be to alter the shape of the interface slightly, such that it would be flatter in 
regions of higher liquid solute concentrations at the interface. 
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CONCLUSIONS 


A series of fully transient simulations of horizontal Bridgman crystal growth under microgravity 
conditions have been performed. The pure bismuth simulation was found to agree reasonably well 
with a similar simulation of the process (Yao et al., 1997), in terms of the convection level and the 
interface curvature. 

For the dilute alloy simulation (Bi-0. 1 at.% Sn), a single dominant counter-clockwise convective 
cell is present for the entire duration of the process. As time proceeds, solute is rejected at the 
interface and the level of solute near the interface increases. Due to the presence of convection, 
radial solute segregation occurs with preferentially higher values at the bottom of the solid. The 
level of segregation increases slightly with time. The maximum values of concentration are small 
and so the assumption of constant alloy melting temperature is realistic. 

For the richer alloy (Bi-1 .0 at.% Sn), the convective field is much more complex. Higher levels 
of solute rejection at the interface cause higher levels of solutal convection. Initially, a single, 
thermally driven, counter-clockwise rotating cell is present. As time proceeds and solute 
accumulates at the interface, a secondary, solute-driven clockwise rotating cell develops. The effect 
of this convective pattern is to yield significant levels of segregation in the solute at the interface 
(and hence in the solid). This segregation is such that the values towards the top of the domain are 
highest. The higher values of concentration near the interface would result in a significant change in 
the melting temperature for the alloy. For this richer alloy, solute build-up is significant and so the 
concentration-dependence of the melting temperature becomes important. 

Finally, simulations of Bridgman crystal growth in microgravity, incorporating the effects of a 
concentration-dependent melting temperature were performed. The results agree well with data 
obtained from the flight experiment. The combination of numerical simulations and experimental 
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data presented here has described the role of thermosolutal convection and solute partitioning during 
directional solidification. 
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Table 1. Thermophysical Properties. 


The properties (except the solutal diffusion coefficient D) for solid and liquid bismuth were 
evaluated at the mean solid and the mean liquid temperatures of 160.7 and 485.7* C, respectively. 
The solutal diffusion coefficient was evaluated at a position halfway into a typical solute boundary 
layer, since solute diffusion occurs mainly in this region. The reference density was considered to be 
equal for both phases. 


Property 

Units 

Value 

Reference density, * o 


10070 

Thermal conductivity of the solid, ks 

W/mK 

6.872 

Thermal conductivity of the liquid, kL 

W/mK 

14.66 

Specific heat of the solid, c p s 

J/kgK 

132.6 

Specific heat of the liquid, c p l 

J/kgK 

135.3 

Enthalpy of freezing, • H 

kJ/kg 

52.3 

Viscosity of the liquid, p 

Ns/m 


Melting temperature, T m 


271.3 

Thermal expansion coefficient, • t 



Thermal conductivity of the (fused silica) ampoule, k w 

W/mK 

2.10 

Density of the (fused silica) ampoule, • w 

mm 

2020 

Specific heat of the (fused silica) ampoule, c pw 

J/kgK 

1066.8 

Partition (or segregation) coefficient, k p 

— 

0.029 

Diflusivity of liquid Sn in liquid Bi, D 

| 

iiiiggymn 

Solutal expansion coefficient, • c 

| 

• 0.305 

Prandtl number, Pr 

— 

0.01144 

Lewis number, Le 

— 

3074 

Slope of liquidus line, m 

K/(vol%) 

-2.547 














































List of Figures 

Fig. 1. Variation of the activation energy for creep, normalized by the activation energy for 
lattice self-diffusion, with the homologous temperature for several materials showing the 
probable dominance of dislocation core and lattice self-diffusion at intermediate and high 
temperatures, respectively [22]. 

Fig. 2. Dependence of the experimental activation energy for creep for (a) & (c) NaCl single 
crystals [40,41], (b) Cu [34] and (d) NiAl [46] on either the homologous temperature or the 
normalized stress. Note that E is the Young’s modulus in (d). The horizontal broken lines 
represent the activation energies for lattice diffusion and pipe diffusion, Q p , of the different 
species. Region I in (d) represents a low stress regime, whereas II and III represent the power-law 
and the exponential creep regimes, respectively. The experimental data are also compared with 
the predictions of pipe-to-lattice self-diffusion transition models [22], cross-slip [27,29] and 
obstacle-controlled [6] creep theories. 

Fig. 3. Normalized creep rate vs. normalized stress plots for (a) f.c.c. metals compiled from the 
literature [15], (b) NaCl single crystals [47] and (b) NiAl [46] showing the transition from 
power-law to exponential creep. 

Fig. 4 Effect of grain size on the compressive deformation behavior of (a) FeAl [49,50] and (b) 

& (c) as-cast and powder-metallurgy processed NiAl [48] showing grain size strengthening in 
intermetallic alloys at high temperatures. 

Fig. 5 Creep substructures observed in the secondary creep regime in (a) pure A1 [52], (b) NaCl 
single crystals [47], (c) NiAl [46] and (d) a Ti 47 Al 51 Mn 2 alloy [69]. 

Fig. 6 A simplified schematic diagram showing how an initially homogeneous distribution of 
dislocations could bifurcate into several possible dislocation substructures during deformation. 

Fig. 7 Schematic plots demonstrating that the transient primary creep stages, A, B, C and D 
shown in (a) correspond to the structure-dependent (i.e. S f ) exponential curves in the normalized 
strain rate-normalized stress plot shown in (b). 

Fig. 8 (a) Variation of creep strain, creep rate and the applied stress with time. A homogeneous 
distribution of dislocations is assumed to form as soon as plastic deformation commences, which 
is assumed to evolve to a mixture of light patches (i.e. recovered regions with a lower dislocation 
density), tangles and homogeneous dislocations as shown in (b) at the resolution limit, (c) 
Schematic showing the concept of “hard” and “soft” regions with respect to the light patches 
assumed by the model. 
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